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ABSTRACT 

The trispectrum of the cosmic microwave background can be used to assess 
the level of non-Gaussianity on cosmological scales. It probes the fourth or- 
der moment, as a function of angular scale, of the probability distribution 
function of fluctuations and has been shown to be sensitive to primordial 
non-gaussianity, secondary anisotropics (such as the Ostriker-Vishniac effect) 
and systematic effects (such as astrophysical foregrounds). In this paper we 
develop a formalism for estimating the trispectrum from high resolution sky 
maps which incorporates the impact of finite sky coverage. This leads to a 
series of operations applied to the data set to minimize the effects of contami- 
nation due to the Gaussian component and correlations between estimates at 
different scales. To illustrate the effect of the estimation process, we apply our 
procedure to the BOOMERanG data set and show that it is consistent with 
Gaussianity. This work presents the first estimation of the CMB trispectrum 
on sub-degree scales. 

Key words: cosmic microwave background- statistics 



1 INTRODUCTION 

The Cosmic Microwave Background (CMB) has become the observational tool of excellence 
for probing the statistical nature of inhomogeneities in the universe. The small deviations 
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from homogeneity which have been detected by over two dozen different experiments can be 
directly related to the primordial origin of perturbations in the early universe and therefore to 
fundamental physics at very high energies. A new threshold was crossed in the experimental 
forum with the high resolution, high sensitivity mapping of significant fractions of the CMB 
sky by the BOOMERanG (de Bernardis et al. 2000) and MAXIMA (Hanany et al. 2000) 
experiments. A careful analysis of the variance of fluctuations in these maps has led to 
accurate estimates of the angular power spectrum, far surpassing previous experimental 
analyses on small angular scales. More recent results from BOOMERanG (Netterfield et al. 
2002, Ruhl et al. 2002), MAXIMA (Lee et al. 2001), and from other experiments like DASI 
(Halverson et al. 2002), CBI (Pearson et al. 2002), VSA (Grainge et al. 2002), ACBAR (Kuo 
et al. 2002) and ARCHEOPS (Benoit et al. 2002) have posed our knowledge of the CMB 
angular power spectrum on even more solid ground. However, there is more information in 
the CMB fluctuations than what is provided by its power spectrum alone. The standard 
way to extract this information is to analyse high signal to noise maps of the CMB field like 

those already produced by BOOMERanG or the ones that the MAP I satellite is expected 
to provide shortly. 

There is, therefore, strong motivation to attempt a more detailed study of the CMB 
sky; in principle one would like a complete characterisation of the probability distribution 
function of the fluctuations in the CMB with the hope that it might probe more fundamental 
features of the origin of structure in the universe. For example one relatively stringent test 
of whether the origin of fluctuations is due to a standard, single-fleld inflationary model 
is if the CMB is a realization of a nearly Gaussian distribution (Coles and Barrow 1987), 
while other models like the curvaton might lead to a much larger non- Gaussian contribution 
(Lyth & Wands 2002; Bartolo & Liddle 2002; Bernardeau & Uzan 2002; Bernardeau & 
Uzan 2002). Many secondary anisotropics hke the Ostriker-Vishniac (OV) effect (Ostriker & 
Vishniac 1986) and the Sunyacv-Zeldovich effect (Sunyacv & Zeldovich 1980) can introduce 
measurable non-Gaussianity while foregrounds and systcmatics may contribute as well. 

A standard method of parametrising an arbitrary probability distribution function is in 
terms of all its higher order moments. In the case of statistically homogeneous and isotropic 
flelds it is more convenient to consider moments of the Fourier transform of the fleld; these 
symmetries will impose a set of selection rules which pick out the true degrees of freedom. The 
past few years have seen initial attempts at constraining these moments by measuring them 
with the available data. A series of analyses of the bispectrum (the cubic moment) of the 
COBE data have revealed the presence of a non-Gaussian systematic (Ferreira, Magueijo, & 
Gorski 1998; Banday, Zaroubi & Gorski 2000; Komatsu ct al. 2001). This has been confirmed 
with an analysis of the trispcctrum, the quartic moment (Kunz et al 2001; Komatsu 2002). 
But no primordial non-Gaussianity was detected. On smaller angular scales analyses of the 
QMAP and QMASK(Park et al. 2001) and MAXIMA (Wu et al 2001; Santos et al 2002) 
data have shown that their observations are consistent with the assumption that the CMB 
anisotropics are the result of an isotropic Gaussian random process. The analyses of these 
data sets have also revealed that statistical methods may be sensitive to the data processing 
pipehne. Moreover, different technical issues must be confronted if one is considering finite 
sky coverage as opposed to full sky coverage. 

The trispectrum, which we consider in this study, probes a rather different kind of non- 
Gaussianity than the bispectrum; (Aghanim et al 2003) found it to be very sensitive to point 
sources and (Castro 2003) showed it to be a powerful probe of the OV effect. It can also be 

t http://map.gsfc.nasa.gov/ 
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used to detect weak lensing in the CMB (Cooray & Kesden 2002). In particular, as described 
in (Bernardeau 1997), weak lensing does not introduce a three-point correlation function, 
meaning that its expected bispectrum is zero. The first non trivial higher order correlation 
function to detect weak lensing effects is the trispectrum. Using the trispectrum to test for 
non-Gaussianity in high-resolution CMB maps complements therefore other analyses using 
the bispectrum, like (Santos et al 2002). The purpose of this paper is to improve, extend 
and apply the method for estimating the trispectrum first proposed in Kunz et al (2001) 
where it was applied to the full sky map produced by the COBE satellite. 

Using techniques developed for performing the operation in pixel space (Ferreira, Maguei- 
jo, & Silk 1997; Spergel & Goldberg 1999; Hu 2001) we present a formalism which can 
easily be extended to high is and discuss a robust method for identifying an orthogonal 
set of estimators in the full-sky case. Wc then discuss the various numerical and statistical 
problems one faces when analysing finite sky coverage. Wc finally use the BOOMERanG 
data as a test case to extract the first estimate of the trispectrum on sub-degree angular 
scales. 



2 FORMALISM AND NOTATION 

In this section we present the notation that will be used throughout this work. We start 
with a temperature anisotropy field defined on the sphere, -^(n); it may be expanded in 
terms of spherical harmonic functions, Y^^(n): 

AT 

— {n)=Y,aimYim{n). (1) 
em 

For any theory of structure formation, the coefficients are a set of random variables; we 
shall restrict ourselves to theories which are statistically homogeneous and isotropic. The 
power spectrum Cg of the temperature anisotropies is then defined by 

(a^mOl'm') — Ce5u>5mm'- (2) 

If we consider the 3-point function of the temperature field, we obtain the bispectrum, 
defined as 

/ \ ^1 ^2 4 \ o fo^ 

{ai^m^ai^m^ai^mi) = I ^2 mg j ^^i^2^3" 

The term (• • •) is a Wigner 3 J symbol, which arises due to the "selection rules" of the 
moments. 

Following the same steps, we can construct the 4-point function and the associated 
trispectrum. We represent the rotationally invariant solution for the trispectrum as in (Hu 
2001): 

{a4rrn,atrn^.a,,„,,a,^^^ ^ S ( -M ) ( m, M ) (4) 

Using the orthogonality properties of the Wigner 3 J symbols and the relation Q — T + G, 
we can invert the equation (4) to obtain the estimator 



2 L \ ( h h L \ ^ (5) 



m\m2mzm4M \ / \ , 



© 0000 RAS, MNRAS 000, 000-000 



4 G. De Troia et al. 



The term G|^^(L) represents the unconnected Gaussian contribution and it is given in 
(Hu 2001) as: 

G\l\l{L) = (-i)ai+'3)^(2/, + i)(2/3 + 1)Ci,Ci,6i,,AjJlo + (2L + l)Q,Ci, x (6) 

The term T^^^^ (L) is the connected part of the angular trispectrum and its expectation 
value is exactly zero for a Gaussian field. This means that the connected part is sensitive 
to the presence of non Gaussianities. The unconnected term is non-zero only for L = or 
^1 = ^2 = ^3 = ^4, but only with full sky coverage. In the case of incomplete sky coverage 
the unconnected terms can contaminate all other modes. We will discuss this situation in 
section (3). 

For the purpose of this work we have not computed the possible trispectrum compo- 
nents. We concentrated only on the simpler case ii = £2 = h = = i-e. the diagonal 
component. Recent papers have investigated in detail the power of the diagonal trispectrum 
in the presence of some non-Gaussian signals mentioned in the introduction. In particular, 
(Aghanim et al 2003) has shown that for simulated point-source maps the diagonal trispec- 
trum is much more powerful than the nearly diagonal estimator + l,i + 2,i + 3), even 
though the latter does not contain a Gaussian contribution. Also, in (Castro 2003) it is 
discussed how the OV effect generates a signature on the diagonal tispectrum which could 
easily be detected on the arcminute scales probed by the PLANCK-t mission. 

Finally it should be noticed that computing all components of the trispectrum is a serious 
computational challenge. Many modes are also correlated due to the limited sky coverage. 
For these reasons we have decided to restrict this analysis to the case of ii — £2 — '^s — — 

We start with the method described in (Kunz et al 2001). We define 

n 
a=0 

where the a arc the components of the trispectrum that we wish to estimate, and r 
is a tensor which we have to determine in order to construct an estimator for T^^a- The 
geometrical considerations stated above, together with the required symmetries with respect 
to the interchange of {ii,mi} pairs suggest 

rrr:,m„ = ,^J~^^ [7^ m2 M )[-M Tu, m, ) + i^eq.permut. (8) 

AlthoTigh the fs satisfy all the correct symmetries, they define an over-complete basis. To 
correct for this deficiency we define an orthonormalised set of tensors 



r"' = > raa-a,i Iq\ 

mxmivnzra^ / j £ mim2fn3m4' v / 

where the matrix is derived from the required property that the r be orthogonal with 
respect to the product given in equation (5) of Kunz et al (2001). The estimator of the 
trispectrum is then given by 

(10) 

mim2m3m4 

t http: / / astro.estec.esa.nl/SA-general/Projects/Planck/ 
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'^mim2m3m4 ^^"^i ^^"12 '^ims (^irrn 

mi a=Q 

e 

0=0 

Note that there are only int(£/3) independent estimators due to the symmetry properties 
of the a£ra- In this paper we will consider the "normalised" trispectrum used in Kunz et al 
(2001) where we divide each estimate of the trispectrum by (Q)^, where C(, = i^^aim!^. 
Its statistical properties are equivalent to the ones of the unnormalised estimator, and it is 
more robust with respect to fluctuations in the power spectrum ( Aghanim et al 2003) . 



3 APPLICATION TO HIGH RESOLUTION MAPS WITH INCOMPLETE 

SKY COVERAGE 

In this paper we will be focusing on a high-resolution map with incomplete sky coverage, 
in particular the BOOMERanG data set. This leads to a set of algorithmic problems which 
did not have to be addressed in Kunz et al (2001). The three problems we wish to highlight 
are: 

Speed: The numerical evaluation of Wigner 3 J coefficients for large values of i becomes 
time consuming and practically unfeasible. Indeed beyond the COBE resolution of a maxi- 
mum i of approximately 25 it is not possible to estimate the T"'^ sufficiently rapidly for a 
robust Monte Carlo assessment of the statistics. 

Gaussian Contamination: The finite sky coverage will induce correlations between the 
estimators with different values of i and a (or a). As a consequence, all estimators may be 
heavily contaminated by the Gaussian (or disconnected) contributions to the maps. 

Correlations: The correlations between modes in the cut sky mean that the T"'^ will be 
even more correlated than in the full sky case. 

We shall now focus on the solutions we propose to these three problems 



3.1 Speed 

We have opted to use the method described in (Hu 2001; Spergel & Goldberg 1999) for 
calculating T"-'^^: we define a new set of sky maps weighted in rings centred around a point 



'2£+ 1 



JdhT{h)Pe{h-q) (11) 

To implement this method we start with the relation (11) and we use the relation (1) to 
express the temperature T as a function of spherical harmonics and the relation: 

Peii^ ■ q) = ^ E y;min)Ya^) (12) 

to express also the Legendre polynomials as a function of spherical harmonics. Combining 
them with equation (11) we obtain: 



/ An 

e€(q) = J ^jqTY E «^m>€m(q) (13) 



© 0000 RAS, MNRAS 000, 000-000 



6 G. De Troia et al. 



The e£ calculation is quite fast because we can use the fast Fourier transform on rings of 
equal latitude (Muciaccia et al 1997). 

We can then rewrite equation (5) in terms of this new set of sky maps (Komatsu 2002): 



2a 



fa;e__}__ fie* fU (■,A\ 

^" "T M=-2a 



where 

,u _ . /2^ + l ( ^ ^ L 

V 



V 47r I 



I dh[ee{n)e,{h)]Yl^{h) (15) 



If we expand the Wigner 3J symbols in terms of spherical harmonics and use the addition 
theorem we obtain: 

f'^-'' = Nil I'^l ^ee{h)ee{h)ee{q)ee{q)P,^{h ■ q) (16) 
where 

We have thus computed a set of T"'^ which we can orthonormalise to get the estimator for 
the trispcctrum T°'^. This method is very fast, especially when estimating the trispectrum at 
high values of i. Note that direct evaluation of the Wigner 3J coefficients, e.g. by recurrence 
relations, would result in an 0{£^) problem, requiring ~ 10^^ operations for £ ~ 1000. 



3.2 Gaussian Contamination 

Kunz et al (2001) found that the purely Gaussian contribution to the trispectrum (the 
disconnected part) corresponds to the a = term. By orthogonalising all other estimators 
with respect to this tensor it is possible to remove the Gaussian contribution exactly on a 
map by map basis. The resulting estimators are only sensitive to non-Gaussian contributions, 
i.e. in the case of Gaussian skies they would have a zero expectation value. Additionally, the 
variance of this estimator is shown to be minimal (Kunz et al 2001). In the cut sky case, 
there are strong cross correlations between components of the trispectrum with different 
values of £ and a. In this case, the orthogonalisation method fails. Given that the Gaussian 
contribution to the trispectrum may be much larger than the non-Gaussian contribution, it 
is essential that we remove it as completely as possible nonetheless. 

To overcome these problems we have chosen to employ the following Monte Carlo scheme: 
we generate an ensemble of maps with the same angular power spectrum, sky coverage and 
noise as the data maps we want to analyse. We estimate the T"'^ from each map and calculate 
the mean of these quantities over the whole ensemble. Let us denote this mean by f^'^. We 
then use this quantity to correct for the Gaussian contamination in the estimate of the 
trispectrum from the data by defining = T"'^ — Tq^. Note that, by so doing, we are 
removing the Gaussian contamination before performing the full-sky orthogonalisation, i.e. 
before multiplying by C 



3.3 Correlations 

The fact that there is only limited sky coverage also implies that there will be correlations 

between values of the trispectrum at different £s. This is shown in section 5 below. To strongly 
suppress the correlations and to end up with a simple covariance matrix, we consider band 
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averaged values of the trispectrum. Since there is no a priori given band width, we study 
the cases of Ai — 40, 50 and 60 and Aa = 10 and 15. The choice Ai ~ 50 is consistent 
with previous analyses of the BOOMERanG power spectrum (see, e.g., (Ruhl et al. 2002)). 
In this way we can check the sensitivity of the results to the chosen band size. 

4 NUMERICAL IMPLEMENTATION AND CONSISTENCY TESTS 

The process wc use is basically the same as in a number of previous analyses (Kogut et al 
1996; Ferreira, Magueijo, & Gorski 1998; Kunz et al 2001; Santos et al 2002). We generate 
an ensemble of Gaussian maps with the same angular power spectrum and noise property 
as estimated from the BOOMERanG data and the same sky coverage. We then apply our 
estimators to the set of maps to obtain a distribution for each estimator in the Gaussian 
case. In particular we characterise the full distribution in terms of the mean values of the 
estimators and the covariancc matrix between them. These quantities are used to define a 
standard multivariate ^ goodness of fit. The estimators of the trispectrum are then 

evaluated from the BOOMERanG data. In section 5 we will discuss their behaviour. The 
goodness of fit of these estimators is compared against its distribution for a new ensemble 
of Gaussian maps. From this comparison we can quantify the confidence with which the 
data can be said to be Gaussian from the point of view of our estimator. It is clear that the 
numerical details of this process must be well understood if we are to believe in our results. 
We focus on the particularities of the analysis in this paper which differ from previous 
analyses. 

It is important to compare the results using this hybrid pixel/harmonic analysis with the 
standard methods which have been used before. We do so by looking at the two lower order 
statistics, i.e. we have calculated the power spectrum Q and the bispectrum using this 
new approach as well as summing up the 3J symbols. The relevant expressions are: 

Ce = ^ J dci\ee{q)f (18) 



and (Spergel et al. 1999) 




We have compared the Ci and Bm using these expressions with the standard results 
obtained using aim and the Wigner 3J symbols, for a maximum i value of 1000. Using 
a set of CMB Gaussian maps with the best fit power spectrum measured by Boomerang 
(Netterfield et al. 2002) and a pixel resolution of 7' we have found that the bispectrum 
obtained with is affected by a pixelisation effect for high values of £ (while the power 
spectrum shows no difference). To check for this, wc have done the same analysis with 
higher resolution (~ 3 ) and have found that the pixelisation effect vanishes. Given that 
we are restricted to the pixelisation level of the data, we can use the comparison of the 
two estimates of the bispectrum to define a maximum £ out to which wc can trust the new 
estimate of the trispectrum. Note that it is computationally intractable to perform such a 
comparison in terms of the trispectrum, although this would be preferable. From Fig. 1 one 
can see that discrepancies arise for £ > 800 and we chose not to estimate the trispectrum 
beyond £ — 700, leaving a conservative margin as we did not test the trispectrum itself. 
Furthermore, we chose not to consider any £ below 100 because the BOOMERanG data are 
not very sensitive to these modes, due to limited sky coverage and data filtering. 

Another novelty in our analysis (as compared to that of Kunz et al. 2001) is the method 
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Figure 1. In this figure we represent the normaUsed bispectrum /| = Bi/C^^'^ calculated for a full sky Gaussian map with a 
pixel size of 7 . The red crosses show /| obtained with the method, the black boxes /| with and the 3J symbols. We 
can see that at ^ ~ 800 there is a quite evident difference between two plots, due to a pixelisation effect. We limit therefore our 
analysis to I < 700. 

for constructing There, a Gram-Schmidt (GS) procedure was used to calculate the 

orthonormal transformation matrix -C"'". Due to the inherent instability of the GS procedure, 
it is not applicable to large matrices, i.e. for large ^. We have opted therefore to use an 
alternative orthonormalisation method. We subtract the a = part and use a Jacobi routine 
to obtain a spectral decomposition (SD) of the remaining matrix. The eigenvectors of the 
non-vanishing eigenvalues form then the transformation matrices. This method is robust 
and, moreover, gives us an unambiguous procedure for ordering the estimators through the 
different eigenvalues. As a strong consistency test we have applied our trispectrum code to 
the COBE data and compared the results with those of (Kunz et al 2001). The results of the 
SD method lie, up to a possible sign change, very close to the original ones (see figure 2). 
In any case, the statistical significance of the results (and the conclusions one can extract) 
are the same as in Kunz et al (2001). We advocate the use of the SD method from now on, 
even in the case of analyses limited to low ^s. 

For our analysis we have used the best four of the six 150 GHz channels of the BOOMERanG 
1998 flight; we naively coadd the data taken at the scan speed of 1 degree per second (1 
dps). We simulate three sets of 1000 Gaussian maps each. In fact, we need three statistically 
independent ensembles of simulated maps for our analysis: one to estimate the Gaussian 
contribution described above, another to estimate our estimator's covariance matrix and a 
third one to perform the actual comparison with real data. To produce the maps we fol- 
low the very same steps used during the real Boomerang data reduction. To generate a 
map we use timestream simulations created with the actual flight pointing and transient 
flagging. The signal component of these time streams is generated from simulated gaussian 
CMB maps, while the noise component is from gaussian realizations of the measured detec- 
tor noise power spectrum. The noise spectrum is determined using the iterative procedure 
described in Ferreira k Jaffe (2000) and Netterfield et al. (2002). 

To reduce the effects of 1// noise on this naively binned map, a brick- wall highpass 
Fourier filter is first applied to the timestream at a frequency of 0.1 Hz. A notch filter is 
also applied between 8 and 9 Hz to eliminate a non-stationary spectral line in the data. This 
only affects angular scales well above C. = 1000 and is therefore irrelevant for our analysis. 

The coaddition of several channels is achieved by averaging the maps (both from the 
data, and from the Monte-Carlos of each channel). Each channel has a slightly different 
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Figure 2. In this figure we reproduce the COBE results for the trispectrum with the GS and SD orthonormaUsation methods 
(see text) and compare them with (Kunz et al. 2001). The black crosses are the (Kunz et al. 2001) results, green squares are 
GS method results and blue triangles are SD method results. 

beam size, which is taken into account in the generation of the simulated maps. We select 
the most central region of the scan by applying an elliptical mask as in (Netterfield et al. 
2002). This corresponds to ~ 1.8% of the full sky. The mask selects a region of approximately 
uniform coverage and high signal to noise and comprises ~ 57000 pixels of size ~ 7' each 
in the Healpix pixelization scheme (Gorski et al. 1998). We refer the reader to Ruhl et al. 
(2002) for a thorough description of a simulation pipeline (based on the MASTER/FASTER 
algorithms described in (Hivon et al. 2000)) very similar to the one used here. 



5 RESULTS AND APPLICATION TO THE BOOMERANG DATA 

To show in detail the method proposed in Sections 2 and 3 we are going to discuss the results 
obtained at each step from both the data and the Monte-Carlo simulations. We start with 
the estimate of the T"'^ without Gaussian corrections. In the top panel of Figure 3 we plot 
T"'^ as a function of i for selected values of a. We can highlight two features. Firstly the 
T°'^ are highly correlated for adjacent values of i due to the finite sky coverage, as expected. 
Secondly, and because of the finite sky coverage, there is a strong contamination from the 
disconnected component of the trispectrum. This is evident in the fact that the values of 
T"'^ scatter about the (C^)^ and that the 95% confidence limits are not centred about zero. 
As one would expect the lower the value of a, the more contaminated the estimate is by the 
disconnected part. As advocated in Section 3, we correct for the contamination due to the 
disconnected component by using a Monte Carlo ensemble (of 1000 realizations) to generate 
a correction. This can be seen as a bias which must be subtracted off all estimates of T"'^ 
with a > 0. In the bottom panel of Figure 3 we plot the "Gaussian corrected" estimate 
of T"'^ with corresponding 95% confidence limits. As expected the estimates now scatter 
around zero, while the confidence limits, although not necessarily symmetric around the i 
axis are much more centred. The remaining asymmetry is merely a manifestation that for 
low as the distribution of the T"'^ is slightly skewed. 

Given that we will be working with normalised estimates of the trispectrum (as in Kunz 
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Figure 3. Top panel: An estimate of the non-orthogonaliscd trispectrum, f "'^ (multiplied by for a = 5, 100, 400 from the 
BOOMERanG data (crosses) and the corresponding 95% confidence limits from the 1000 Monte-Carlo simulations. Bottom 
panel: An estimate of the non-orthogonalised trispectrum corrected for Gaussian contamination, T"'^ (multiplied by £'^) for 
a = 5, 100, 400 from the BOOMERanG data (crosses) and the corresponding 95% confidence limits from the 1000 Monte-Carlo 
simulations. 




Figure 4. An estimate of the non-orthogonalised trispectrum corrected for Gaussian contamination and normalised f "'^/(C^)'^ 
for a = 5, 100,400 from the BOOMERanG data (crosses) and the corresponding 95% confidence limits from the 1000 Monte- 
Carlo simulations. 

et al 2001) it is illustrative to plot the £ dependence of f"''^/ (Q)^ for a few values of a. We 
do this in Figure 4. One can see the dependence on a of the 95% confidence intervals, i.e. 
the i value of minimum scatter depends on a. 

Let us now proceed to the orthornomalised estimators, Ta-f, a selection of estimators are 
plotted for a choice of as in Figure 5. As noted above, the a are limited to a < int(//3), and 
we see a clear suppression of the high a values for each £ (or, in the case of the figure, of 
the low i values for fixed a) , as the maps with limited sky coverage contain less information 
than full sky maps. 
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Figure 5. Top panel: An estimate of the orthogonalised trispectrum corrected for Gaussian contamination and normalised 
Ta;e/{Ci)^ for a = 5,100,200 from the BOOMERanG data (crosses) and the corresponding 95% confidence limits from the 
1000 Monte-Carlo simulations. Bottom panel: The absolute value of the estimate of the orthogonalised trispectrum corrected 
for Gaussian contamination and normalised T„,.^i,j{C()'^ for a = 5,100,200 from the BOOMERanG data (crosses) and the 
corresponding 95% confidence limits from the 1000 Monte-Carlo simulations. 
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Figure 6. The distribution of Monte-Carlo simulated maps (histogram) and data value (vertical line) for the trispectrum 
estimator \Ta-^i\ in the case of Aa = 10 and A/ = 40, 50, 60 (top) and in the case of Aa = 15 and Ai = 40, 50, 60 (bottom) 
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Figure 7. The distribution of Monte-Caxlo simulated maps (histogram) and data value (vertical line) for the trispectrum 
estimator T^.^ in the case of Ao = 10 and Al = 40, 50, 60 (top) and in the case of Ao = 15 and Al = 40, 50, 60 (bottom) 
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Table 1. Probability that the Gaussian models have a greater than the data value for both the trispectrum estimators and 
for different bin widths in i and a. 

Once we have calculated the trispectrum estimators both for the BOOMERanG data and 
for the Monte-Carlo simulations, we can proceed to get the distribution for the simulated 
Gaussian maps and compare it to the data. We used two different approaches: one taking 
as estimator Ta-/ (the orthogonalised trispectrum corrected for Gaussian contamination and 
normalised to Cf) and the other one taking its absolute value \Ta-e\- 

We construct a standard multivariate 

= E i{Tl;a)G - %a)Cll,y{{f,,,)G ' f,.^,) (20) 

t,e',a,a' 

deriving the expectation values {Ti-a)G and the covariance matrix C^^^'^' = {T^aTf .a_')G ~ 
^i;a} g{Ti' -a') G ffom One of the two remaining Monte-Carlo ensembles. As discussed earlier, 
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we don't sum over all i and a, but bin both a and i, varying the size of the bins. Finally we 
calculate the distribution from the last ensemble of Gaussian maps. 

In figures 6 and 7 wc show the distribution derived from 1000 Gaussian realizations 
compared to the BOOMERanG data for both estimators and for different bin widths. The 
probability P{x^ > Xb) ^^^^ ^ Gaussian map has a larger than the Boomerang map 
is given in table 1. Although the values vary considerable with the choice of bin- widths, 
none of them is below 5% or 2a. We conclude that the trispectrum does not detect any 
non-Gaussianity in the coadded Boomerang 150 GHz maps. 

6 CONCLUSIONS 

We have applied an improved version of the method of (Kunz et al 2001) for measuring the 

trispectrum to the four best 150 GHz Boomerang maps. To this end, we used maps containing 
only one multipole each to avoid computing the Wigner 3J symbols and subtracted the 
average Gaussian contribution using an ensemble of simulated maps. We then orthogonalised 
these maps and normalised them to We then binned the resulting trispectrum values 
with a variety of different bin sizes, and computed the value, using a full covariance matrix 
estimated from a second ensemble of simulated Gaussian maps. When comparing the data 

value to the Gaussian realizations (obtained from a third ensemble of simulated Gaussian 
maps) we concluded that the trispectrum does not detect any deviations from Gaussianity. 

This work complements the pixel-space analysis (Polenta et al. 2002) and the bispectrum 
analysis (Contaldi et al. in preparation) of the BOOMERanG data. 

In this paper we have studied for the first time the trispectrum of real CMB data with 
sub-degree resolution. The main problem encountered was the limited sky coverage, which 
introduces strong correlations, and prevents the use of orthogonalisation to remove the 
Gaussian (connected) part of the trispectrum. We expect therefore that the MAP and Planck 
satellites will be able to improve on this analysis considerably, but this study provides a proof 
of feasibility for measuring the trispectrum of full sky high resolution maps as well as first 
results on small angular scales. 
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